function ObjectSegmentation(nfindobjs, infile, fin, fout, fbackground, fsubtract, handles)


flen = length(infile);
f = infile(1:flen-4);


templatefile = strcat(f, '_templates.mat');

load(fin);

v = mmreader(infile);
totalframes = v.NumberOfFrames;


load(fsubtract);
cc = [cr1, cr2, cg1, cg2, cb1, cb2, ch1, ch2, cs1, cs2, cv1, cv2];

load(fbackground)

objs1 = find(objs == 1);

bCreateFile = 1;

if bCreateFile == 1
    bsubtract = 0;
    outfile = templatefile;
    ntrainimgs = 30;
    maxbestscore = 90;
    minimgs = 5;
    nmin = 20;
    
    n = CreateTrainingSet(v, fin, fbackground, bsubtract, outfile, ntrainimgs, nfindobjs, maxbestscore);
        
    while n < minimgs & maxbestscore > nmin
        
        maxbestscore = maxbestscore-nmin;
        n = CreateTrainingSet(v, fin, fbackground, bsubtract, outfile, ntrainimgs, nfindobjs, maxbestscore);
        
    end            
end

load(templatefile)

a = xlengths;
b = ylengths;


xsort = sortrows([a' b'], 1);
ysort = sortrows([a' b'], 2);

p = 0.05;
len = length(a);
n = floor(len*p+1);


if len == 1
    return
end

xbottratioavg = mean(xsort(1:n,1)./xsort(1:n,2));
xbottratiostd = std(xsort(1:n,1)./xsort(1:n,2));

xtopratioavg = mean(xsort(len-n:len,1)./xsort(len-n:len,2));
xtopratiostd = std(xsort(len-n:len,1)./xsort(len-n:len,2));

ybottratioavg = mean(ysort(1:n,2)./ysort(1:n,1));
ybottratiostd = std(ysort(1:n,2)./ysort(1:n,1));

ytopratioavg = mean(ysort(len-n:len,2)./ysort(len-n:len,1));
ytopratiostd = std(ysort(len-n:len,2)./ysort(len-n:len,1));


xbottmax = xsort(n,1);
ybottmax = ysort(n,2);
xtopmin = xsort(len-n,1);
ytopmin = ysort(len-n,2);

c = 1.0;
xbottratiomin = xbottratioavg-(c*xbottratiostd);
ybottratiomin = ybottratioavg-(c*ybottratiostd);


xtopratiomin = xtopratioavg-(c*xtopratiostd);
ytopratiomin = ytopratioavg-(c*ytopratiostd);


len = length(objs1);



btemplate = 0;


%Identification and template matching constants
mincorr = 0.25;
nstep = 5;

pmin = 0.2;
pmax = 0.9;

pminlargearea = 0.65;
pminpixelslarge = 0.75;

pminpixels = 0.55;
pminxyaxis = 0.8; %0.65;

pminpixelsoftotal = 0.8;

pmaxoverlap = 0.4;
pmaxratio = 4.0;

pmintemplatematchpixels = 0.6;

ncount = 1;

nmaxsametemplate = 60;
n = 0;

nMaxRuns = 1;

global bRunCcode;

bRunCcode = 1; %bRunCcode = 1 -> run c-code, bRunCcode = 0 -> run matlab code


b = find(objs == nfindobjs);

eas = veas(b,:);
ebs = vebs(b,:);

meaneas = mean(eas(:));
stdeas = std(eas(:));

meanebs = mean(ebs(:));
stdebs = std(ebs(:));

eabs = eas(:).*ebs(:);
meaneabs = mean(eabs);
stdeabs = std(eabs);

c = 2.0;
c1 = 2.0;
c2 = 0.5;

minea = meaneas - stdeas*c;
maxea = meaneas + stdeas*c;
mineb = meanebs - stdebs*c;
maxeb = meanebs + stdebs*c;

minea1 = meaneas - stdeas*c2;
maxea1 = meaneas + stdeas*c2;
mineb1 = meanebs - stdebs*c1;
maxeb1 = meanebs + stdebs*c1;

minea2 = meaneas - stdeas*c1;
maxea2 = meaneas + stdeas*c1;
mineb2 = meanebs - stdebs*c2;
maxeb2 = meanebs + stdebs*c2;

mineos = 0.15;
minpscore = 0.1;

cthresh = 90;

bsegmented = zeros(len, 1);

%starting frame for segmentation, don't segment until all objects have 
%first been found as post-processing will not have been applied until then
sf = find(objs == nfindobjs, 1, 'first'); 
startframe = sf(1) + 1;

stepcount = 1;

%Go through array of frames where less than the number of objects searched
%for was found (e.g looking for 2 objects but only 1 was found). In this
%case, check if each identified object is actually 2 or more objects.  
tic
for i = startframe:len %45ii %1:len %1:130 %1:400 %1:len        
    

    templateobjs = 2;
        
    index = objs1(i);    
    vindex = index + nprocessstart-1;   
    
    [i, vindex]
            
    img = read(v,vindex);

    vx1 = vx1s(index, 1); 
    vx2 = vx2s(index, 1); 
    vy1 = vy1s(index, 1); 
    vy2 = vy2s(index, 1); 
        
    if vx1 > 0 & vx2 > 0 & vy1 > 0 & vy2 > 0
        mainimg = img(vx1:vx2, vy1:vy2, :);                             

        %Get binary image of current frame (which has less than the number of
        %objects we're looking for)
        mainbwimg = GetBWArea(img, fbackground, vx1, vx2, vy1, vy2, cc, nbwarea, nstremove, nstfill, dthresh);         
        tempimg = mainimg;
        tempbw = mainbwimg;
           
        %If a template match was performed with the previous frame(s) then 
        %check if the correlation between the current frame and previous
        %frame. If the previous and current frame are similar enough then
        %use the previous template again instead of searcing through all 
        %templates to find suitable match
        c = 1000;
        if btemplate == 1            
            msz = size(mainimg);                          
            pm = imresize(previmg, [msz(1) msz(2)]);   
            
            [bestx1, bestx2, besty1, besty2, bestscore] = TempMatchSameSizeImages(mainimg, pm); %template match images of same size
            
            c = bestscore;            
        end    
    
        previmg = mainimg;    
    
        bMatch = 0;
    
        if c < cthresh 
                        
            [o1x1, o1x2, o1y1, o1y2, bestscore1] = TempMatch(mainimg, templateimg1); %scan through image with template and perform template matching
       
            if o1x1 > 0
                mainimg(o1x1:o1x2, o1y1:o1y2, :) = 255;                   
                
                [o2x1, o2x2, o2y1, o2y2, bestscore2] = TempMatch(mainimg, templateimg2);
            
                if o2x1 > 0
                    mainimg(o2x1:o2x2, o2y1:o2y2, :) = 255; 
                                                
                    [o1x1, o1x2, o1y1, o1y2, o2x1, o2x2, o2y1, o2y2, paxis, bMatch, bRunAgain] = FitAndCheckTemplate(o1x1, o1x2, o1y1, o1y2, o2y1, o2y2,...
                    o2x1, o2x2, vx1, vx2, vy1, vy2, tempbw, xbottmax, ybottmax, xtopmin, ytopmin, xbottratiomin, ybottratiomin,... 
                    xtopratiomin, ytopratiomin, pmaxoverlap, pminpixels, pminlargearea, pminpixelslarge, pmaxratio, pminxyaxis, pminpixelsoftotal);    
                                
                    n = n + 1;
                                
                    timg1 = tempbw(o1x1:o1x2, o1y1:o1y2);
                    timg2 = tempbw(o2x1:o2x2, o2y1:o2y2);                                
                 
                    [ex0o1, ey0o1, ephio1, eao1, ebo1, eso1, eto1, eobjs1] = GetEllipse(timg1, o1x1, o1y1);
                    [ex0o2, ey0o2, ephio2, eao2, ebo2, eso2, eto2, eobjs2] = GetEllipse(timg2, o2x1, o2y1);                 
                                    
                    if (eao1 < minea1 & ebo1 > maxeb1) | (eao2 < minea1 & ebo2 > maxeb1)
                        bMatch = 0;
                    end
   
                    if (ebo1 < mineb2 & eao1 > maxea2) | (ebo2 < mineb2 & eao2 > maxea2)
                        bMatch = 0;
                    end                                
                else
                    bMatch = 0;
                end
                
            else
                bMatch = 0;
            end           
        end
        

        if bMatch == 0
        
            bMatch = 0;
            nRun = 0;
       
            ns = 0;
            skipindexes = [];       
              
            while bMatch == 0 & nRun < nMaxRuns     
           
                                        
                o2x1 = 0;           
                templateobjs = 2;             
                          
				%Takes each template from template database and scans through image
                [o1x1, o1x2, o1y1, o1y2, bestscore1, bestindex1] = TempMatchDBImagesWithImage(mainimg, fullimgs);
                
                if o1x1 > 0                                
                    templateimg1 = cell2mat(fullimgs(bestindex1));       
                    mainimg(o1x1:o1x2, o1y1:o1y2, :) = 255;       
                    mainbwimg(o1x1:o1x2, o1y1:o1y2) = 0;                                                              
                else
                    templateobjs = 0;
                end      
            

                if templateobjs == 2                            
                                                        
                    %Takes each template from template database and scans through image but also checks if a new match overlaps with previous match                                                    
                    [o2x1, o2x2, o2y1, o2y2, bestscore2, bestindex2] = TempMatchDBImagesWithImageCheckOverlap(mainimg, fullimgs, o1y1, o1y2, o1x1, o1x2, pmaxoverlap,i);
                       
                    if o2x1 > 0
                        templateimg2 = cell2mat(fullimgs(bestindex2));                    
                    else
                        templateobjs = 0;
                    end
                end          
                
      
                if templateobjs == 2     
                                                                            
                    [o1x1, o1x2, o1y1, o1y2, o2x1, o2x2, o2y1, o2y2, paxis, bMatch, bRunAgain] = FitAndCheckTemplate(o1x1, o1x2, o1y1, o1y2, o2y1, o2y2,...
                    o2x1, o2x2, vx1, vx2, vy1, vy2, tempbw, xbottmax, ybottmax, xtopmin, ytopmin, xbottratiomin, ybottratiomin,... 
                    xtopratiomin, ytopratiomin, pmaxoverlap, pminpixels, pminlargearea, pminpixelslarge, pmaxratio, pminxyaxis, pminpixelsoftotal);            
                                                                                        
                    timg1 = tempbw(o1x1:o1x2, o1y1:o1y2);
                    timg2 = tempbw(o2x1:o2x2, o2y1:o2y2);
                            
                    [ex0o1, ey0o1, ephio1, eao1, ebo1, eso1, eto1, eobjs1] = GetEllipse(timg1, o1x1, o1y1);
                    [ex0o2, ey0o2, ephio2, eao2, ebo2, eso2, eto2, eobjs2] = GetEllipse(timg2, o2x1, o2y1);
                         
                    if ex0o1 == 0 | ex0o2 == 0
                        bMatch = 0;
                    end
                    
         
                    if (eao1 < minea1 & ebo1 > maxeb1) | (eao2 < minea1 & ebo2 > maxeb1)
                     bMatch = 0;
                    end
   

                    if (ebo1 < mineb2 & eao1 > maxea2) | (ebo2 < mineb2 & eao2 > maxea2)
                        bMatch = 0;
                    end  

                else 
                    bMatch = 0;                
                end                                
            
                if bMatch == 0
                
                    if o1x1 > 0
                        ns = ns+1;
                        skipindexes(ns) = bestindex1;
                    end
                    
                    if o2x1 > 0
                        ns= ns+1;
                        skipindexes(ns) = bestindex2;
                    end
                                                            
                    mainimg = tempimg;
                    mainbwimg = tempbw;
                end            
            
                nRun = nRun + 1;
            end
       
        end
       
     
        btemplate = 0;
    
        if bMatch == 1         
                   
            if nRun == 1 %& paxis > pminxyaxis
                btemplate = 1;
            end            
        
            objs(index) = 2;    
                
            vx1s(index, 1) = o1x1+vx1; 
            vx2s(index, 1) = o1x2+vx1;  
            vy1s(index, 1) = o1y1+vy1; 
            vy2s(index, 1) = o1y2+vy1;
        
            vx1s(index, 2) = o2x1+vx1; 
            vx2s(index, 2) = o2x2+vx1;  
            vy1s(index, 2) = o2y1+vy1; 
            vy2s(index, 2) = o2y2+vy1;
    
            vex0s(index,1) = ex0o1+vy1;
            vey0s(index,1) = ey0o1+vx1;
            vephis(index,1) = ephio1;
            vets(index,1) = eto1;
            veas(index,1) = eao1;
            vebs(index,1) = ebo1;
                
            vex0s(index,2) = ex0o2+vy1;
            vey0s(index,2) = ey0o2+vx1;
            vephis(index,2) = ephio2;
            vets(index,2) = eto2;
            veas(index,2) = eao2;
            vebs(index,2) = ebo2;
            
            
            bsegmented(index) = 1;            
                    
         elseif templateobjs > 0              
            ncount = ncount + 1;    
        end         
                  
    
        if btemplate == 0
            n = 0;
        elseif n > nmaxsametemplate
            btemplate = 0;
            n = 0;
        end   
    
    end
        
    stepcount = stepcount + 1;
    p = floor(100*(stepcount/(len+1)));
    s = strcat('Object Segmentation...', num2str(p));
    s = strcat(s, '%');    
    set(handles.textoutput, 'String', s);     
end

toc

save(fout, 'objs', 'vx1s', 'vx2s', 'vy1s', 'vy2s', 'vex0s', 'vey0s', 'vephis', 'veas', 'vebs', 'vets', ...
    'vxheights', 'vyheights', 'vareas', 'nprocessstart', 'imgx', 'imgy', 'imghalfx', 'imghalfy', ...
    'nfindobjs', 'nfindmaxobjs', 'nstartframe', 'nendframe', 'objs1', 'bsegmented');



function [ex0, ey0, ephi, ea, eb, es, et, eobjs] = GetEllipse(img, x1, y1)

CC = bwconncomp(img);
L = labelmatrix(CC);

s  = regionprops(L, 'centroid', 'Orientation', 'MinorAxisLength', 'MajorAxisLength', 'Solidity');
    

cd = cat(1, s.Centroid);
ori = cat(1,s.Orientation);
mil = cat(1,s.MinorAxisLength);
mal = cat(1,s.MajorAxisLength);           
sd = cat(1,s.Solidity);           

index = 1;

if length(sd) > 1
    [m, index] = max(mal.*mil);
elseif length(sd) == 0
    ex0 = 0;
    ey0 = 0;
    ephi = 0;
    ea = 0; 
    eb = 0; 
    es = 0; 
    et = 0; 
    eobjs = 0;
    return
end



ex0 = cd(index, 1)+y1;
ey0 = cd(index, 2)+x1;            
t = ori(index);
ephi = (abs(t-90)/180)*pi;   
et = t;
ea = mil(index)/2;
eb = mal(index)/2;   
es = sd(index);

eobjs = length(sd);
        


function [o1x1, o1x2, o1y1, o1y2, o2x1, o2x2, o2y1, o2y2, paxis, bMatch, bRunAgain] = FitAndCheckTemplate(o1x1, o1x2, o1y1, o1y2, o2y1, o2y2,...
          o2x1, o2x2, vx1, vx2, vy1, vy2, tempbw, xbottmax, ybottmax, xtopmin, ytopmin, xbottratiomin, ybottratiomin,... 
          xtopratiomin, ytopratiomin, pmaxoverlap, pminpixels, pminlargearea, pminpixelslarge, pmaxratio, pminxyaxis, pminpixelsoftotal)


[bint, l, r, t, b, poverlap] = GetOverlap(o1y1, o1y2, o1x1, o1x2, o2y1, o2y2, o2x1, o2x2);

bMatch = 1;
bRunAgain = 1;
paxis = 0;



    

if poverlap < pmaxoverlap 
    
        [o1x1b, o1x2b, o1y1b, o1y2b] = FitSquaresToTemplateObjects(tempbw, o1x1, o1x2, o1y1, o1y2, o2x1, o2x2, o2y1, o2y2);
        [o2x1b, o2x2b, o2y1b, o2y2b] = FitSquaresToTemplateObjects(tempbw, o2x1, o2x2, o2y1, o2y2, o1x1, o1x2, o1y1, o1y2);
        
        o1x1 = o1x1b;
        o1x2 = o1x2b;
        o1y1 = o1y1b;
        o1y2 = o1y2b;
           
        o2x1 = o2x1b;
        o2x2 = o2x2b;
        o2y1 = o2y1b;
        o2y2 = o2y2b;         
        
           
        %Check if objects have expanded such that one object is 'within'
        %the other. 
        if (o1y1 >= o2y1 & o1y2 <= o2y2 & o1x1 >= o2x1 & o1x2 <= o2x2) | (o2y1 >= o1y1 & o2y2 <= o1y2 & o2x1 >= o1x1 & o2x2 <= o1x2) 
            bMatch = 0;
        end
        
    
    
        [o1x1, o1x2, o1y1, o1y2] = CropSquare(tempbw, o1x1, o1x2, o1y1, o1y2);
        [o2x1, o2x2, o2y1, o2y2] = CropSquare(tempbw, o2x1, o2x2, o2y1, o2y2);  
    
        area1 = (o1x2-o1x1)*(o1y2-o1y1);
        area2 = (o2x2-o2x1)*(o2y2-o2y1);
        totalarea = (vx2-vx1)*(vy2-vy1);                               
        
        pixels1 = sum(sum(tempbw(o1x1:o1x2, o1y1:o1y2)));
        pixels2 = sum(sum(tempbw(o2x1:o2x2, o2y1:o2y2)));
        
        totalpixels = sum(sum(tempbw));
        
        ppixels = (pixels1 + pixels2)/ totalpixels;                        
        
        if ppixels < pminpixelsoftotal
            bMatch = 0;
        end     
        
           
        ppixelso1 = pixels1/area1;
        ppixelso2 = pixels2/area2;
    
        pareao1 = area1/totalarea;
        pareao2 = area2/totalarea;                 

        %For larger areas the minimum pixel percentage of the match is more
        %strict since there is a higher chance that it overlaps the other
        %animal (in which case this combination of 2 objects tends to produce
        %a square with larger 'empty areas'). 
        if (pareao1 > pminlargearea & ppixelso1 < pminpixelslarge) | ...
           (pareao2 > pminlargearea & ppixelso2 < pminpixelslarge)    
             bRunAgain = 0;
        end              
        
        
        x1 = min(o1x1, o2x1);
        x2 = max(o1x2, o2x2);
        y1 = min(o1y1, o2y1);
        y2 = max(o1y2, o2y2);

        
        xlen = x2-x1;
        ylen = y2-y1;
        
        vxlen = vx2-vx1;
        vylen = vy2-vy1;
        
        paxis = min(xlen/vxlen, ylen/vylen);                 
else
    bMatch = 0;
end




function [x1,x2,y1,y2] = CropSquare(v, xstart, xend, ystart, yend)

x1 = 10000;
x2 = 0;
y1 = 10000;
y2 = 0;


for i=xstart:xend-2
    for j = ystart:yend-2
        %if v(i,j) > 0  & v(i,j+1) > 0  & v(i+1,j) > 0 & v(i+1,j+1) > 0 %& v(i,j+1) > 0 & v(i,j+2) > 0 & v(i+1,j) > 0 & v(i+1,j+1) > 0 & v(i+1,j+2) > 0 & v(i+2,j) > 0 & v(i+2,j+1) > 0 & v(i+2,j+2) > 0
        if v(i,j) > 0  & v(i,j+1) > 0  & v(i+1,j) > 0 & v(i+1,j+1) > 0 & v(i,j+1) > 0 & v(i,j+2) > 0 & v(i+1,j) > 0 & v(i+1,j+1) > 0 & v(i+1,j+2) > 0 & v(i+2,j) > 0 & v(i+2,j+1) > 0 & v(i+2,j+2) > 0
            
            if i < x1
                x1 = i;
            end
            
            if i+2 > x2
                x2 = i+2;
            end
            
            if j < y1
                y1 = j;
            end
            
            if j+2 > y2
                y2 = j+2;
            end                                    
        end
    end    
end

    
if x1 == 10000 | x2 == 0 | y1 == 10000 | y2 == 0

    x1 = xstart;
    x2 = xend;
    
    y1 = ystart;
    y2 = yend;
end




function [o1x1, o1x2, o1y1, o1y2] = FitSquaresToTemplateObjects(v, o1x1, o1x2, o1y1, o1y2, o2x1, o2x2, o2y1, o2y2)

sz = size(v);
vxlen = sz(1);
vylen = sz(2);

stepstart = 1;
stepsize = 1;
stepend = 1;

xstepmax = vxlen;
ystepmax = vylen;

xlen = xstepmax;
ylen = ystepmax;

bx1 = 1;
bx2 = 1;
by1 = 1;
by2 = 1;


x1 = o1x1;
x2 = o1x2;
y1 = o1y1;
y2 = o1y2;

x1b = o2x1;
x2b = o2x2;
y1b = o2y1;
y2b = o2y2;

x0 = (o1x2 - o1x1)/2.0 + o1x1;
y0 = (o1y2 - o1y1)/2.0 + o1y1;

x0b = (o2x2 - o2x1)/2.0 + o2x1;
y0b = (o2y2 - o2y1)/2.0 + o2y1;

p = 0.75;
xp = (o1x2 - o1x1)*p + o1x1;
yp = (o1y2 - o1y1)*p + o1y1;

xpb = (o2x2 - o2x1)*p + o2x1;
ypb = (o2y2 - o2y1)*p + o2y1;


%Don't expand in a direction if the other object is there
if y1 > y0b
    by1 = 0;
end

if y2 < y0b
    by2 = 0;
end

if x1 > x0b
    bx1 = 0;
end

if x2 < x0b
    bx2 = 0;
end

%If there is overlap between the objects, only
%expand if there is a given amount of overlap 
%(set by p percent overlap above)

%Lower left corner overlaps 
if y1 >= y1b & y1 <= y2b & x2 >= x1b & x2 <= x2b    
    %If object overlaps across point p (for example the middle) of other object 
    %in x or y direction then it is ok to expand in this direction. 
    
    if x2 < xpb %if true, does not overlap over point xp
        bx2 = 0;
    end
    
    if y1 > ypb %if true, does not overlap over point yp
        by1 = 0;
    end
end

%Upper Left corner overlaps
if y1 >= y1b & y1 <= y2b & x1 >= x1b & x1 <= x2b        
    if x1 > xpb 
        bx1 = 0;
    end
    
    if y1 > ypb 
        by1 = 0;
    end
end

%Lower right corner overlaps
if y2 >= y1b & y2 <= y2b & x2 >= x1b & x2 <= x2b
    if x2 < xpb 
        bx2 = 0;
    end
    
    if y2 < ypb 
        by2 = 0;
    end
end

%Upper right corner overlaps
if y2 >= y1b & y2 <= y2b & x1 >= x1b & x1 <= x2b    
    if x1 > xpb 
        bx1 = 0;
    end
    
    if y2 < ypb 
        by2 = 0;
    end
end


for step = stepstart:stepsize:stepend
    
  thresh = step;
  
  %%%%%%%%%%%%%%%%% Object 1
  i = o1x1-step;
  
  if i < 1
      i = 1;
  end
    
  while bx1 == 1 & i >= 1
      
      pixels = sum(sum(v(i:o1x1,o1y1:o1y2)));             
                 
      if pixels >= thresh
         
          o1x1 = i;
      end
      
      i = i - step;
      
  end
  
    
  
  i = o1x2+step;
  
  if i > xlen
      i = xlen;
  end
  
  
    
  while bx2 == 1 & i <= xlen 
      
      pixels = sum(sum(v(o1x2:i,o1y1:o1y2))); 
          
      if pixels >= thresh
          o1x2 = i;
      end
      
      i = i + step;
      
  end
     
  
  i = o1y1-step;
  
  if i < 1
      i = 1;
  end
  

    
  while by1 == 1 & i >= 1 
            
      pixels = sum(sum(v(o1x1:o1x2,i:o1y1)));     
      
      if pixels >= thresh
          o1y1 = i;
      end
      
      i = i - step;
  end
    
  
  i = o1y2+step;
  
  if i > ylen
      i = ylen;
  end
    
  
  while by2 == 1 & i <= ylen 
      
      pixels = sum(sum(v(o1x1:o1x2,o1y2:i)));   
        
      if pixels >= thresh
          o1y2 = i;
      end
      
      i = i + step;
      
  end
      

end






function bw = GetBWArea(img, fbackground, x1, x2, y1, y2, cc, nbwarea, nstremove, nstfill, dthresh)


load(fbackground)
img0 = backgroundimg; 

img0 = img0(x1:x2, y1:y2,:);
img = img(x1:x2, y1:y2,:);

himg0 = rgb2hsv(img0); 
dhimg0 = double(himg0);
dimg0 = double(img0);

[bw, idif] = ProcessImage(img, dimg0, dhimg0, cc, nstremove, nstfill, nbwarea, dthresh);









function [bestx1, bestx2, besty1, besty2, bestscore] = TempMatchSameSizeImages(objectimg, templateimg)

dobjectimg = double(objectimg);

bestscore = 100000;
bestx1 = 0;
bestx2 = 0;
besty1 = 0;
besty2 = 0;



sz = size(templateimg);
templatearea = sz(1)*sz(2);       
                    
dtemplateimg = double(templateimg);
                            
%TempMatchSameSizeImagesC performs template matching between 2 images that are the same size,
%while TempMatchC scans through an image with a template (but the matching 
%itself is performed that same in the 2 functions)                     
[x1, x2, y1, y2, score, nremoved] = TempMatchSameSizeImagesC(dtemplateimg, dobjectimg);
                        
pscore = score/templatearea;
                       
if pscore < bestscore   
    
    bestscore = pscore;        
            
    bestx1 = x1;
    bestx2 = x2;
    besty1 = y1;
    besty2 = y2;
        
end    




function [bestx1, bestx2, besty1, besty2, bestscore] = TempMatch(objectimg, templateimg)

dobjectimg = double(objectimg);

bestscore = 100000;
bestx1 = 0;
bestx2 = 0;
besty1 = 0;
besty2 = 0;

minsize = 0.5;

objsize = size(objectimg);

objarea = objsize(1)*objsize(2);

sz = size(templateimg);
templatearea = sz(1)*sz(2);
        
psize = templatearea/objarea;
        
        
if psize < minsize
                    
    dtemplateimg = double(templateimg);
                            
     %TempMatchC scans through an image with a template while TempMatchSameSizeImagesC performs 
     %template matching between 2 images that are the same size, (but the matching itself is performed the same in the 2 functions)                     
    [x1, x2, y1, y2, score, nremoved] = TempMatchC(dtemplateimg, dobjectimg);
                        
    pscore = score/templatearea;
                       
    if pscore < bestscore   
    
        bestscore = pscore;        
            
        bestx1 = x1;
        bestx2 = x2;
        besty1 = y1;
        besty2 = y2;
        
    end    
end





function [bestx1, bestx2, besty1, besty2, bestscore, bestindex] = TempMatchDBImagesWithImageCheckOverlap(objectimg, fullimgs, o1y1, o1y2, o1x1, o1x2, pmaxoverlap,k)

dobjectimg = double(objectimg);

bestindex = 0;
bestx1 = 0;
bestx2 = 0;
besty1 = 0;
besty2 = 0;
bestscore = 10000000;

objsize = size(objectimg);

len = length(fullimgs);


for i = 1:len  %number of templates                   
        
        templateimg = cell2mat(fullimgs(i)); 
                        
        tempsize = size(templateimg);
                
        bCompare = 1;
        
        %Resize if needed, and set bCompare = 0 if template image is bigger
        %than object image (which contains 2 objects while templte image
        %has only 1).
        if tempsize(1) >= objsize(1) & tempsize(2) < objsize(2)         
            templateimg = imresize(templateimg, [objsize(1)-1 tempsize(2)]);            
        elseif tempsize(1) < objsize(1) & tempsize(2) >= objsize(2)
            templateimg = imresize(templateimg, [tempsize(1) objsize(2)-1]);
        elseif tempsize(1) >= objsize(1) & tempsize(2) >= objsize(2)
            bCompare = 0;
        end
 
        if bCompare == 1
        
            dtemplateimg = double(templateimg);

            [x1, x2, y1, y2, score, nremoved] = TempMatchC(dtemplateimg, dobjectimg);
            
            %template size may have changed so get again
            tempsize = size(templateimg);
            templatearea = tempsize(1)*tempsize(2);
        
            pscore = score/templatearea;        
            
            

            if pscore < bestscore & x1 > 0 
                
                [bint, l, r, t, b, poverlap] = GetOverlap(o1y1, o1y2, o1x1, o1x2, y1, y2, x1, x2);                             
                
                if poverlap < pmaxoverlap
                
                    bestscore = pscore;
                    bestindex = i;
            
                    bestx1 = x1;
                    bestx2 = x2;
                    besty1 = y1;
                    besty2 = y2;          
                end
            end   
        end
end



function [bestx1, bestx2, besty1, besty2, bestscore, bestindex] = TempMatchDBImagesWithImage(objectimg, fullimgs)

dobjectimg = double(objectimg);

bestindex = 0;
bestx1 = 0;
bestx2 = 0;
besty1 = 0;
besty2 = 0;
bestscore = 10000000;

objsize = size(objectimg);

len = length(fullimgs);


for i = 1:len  %number of templates                   
        
        templateimg = cell2mat(fullimgs(i)); 
                        
        tempsize = size(templateimg);
                
        bCompare = 1;
        
        %Resize if needed, and set bCompare = 0 if template image is bigger
        %than object image (which contains 2 objects while templte image
        %has only 1).
        if tempsize(1) >= objsize(1) & tempsize(2) < objsize(2)         
            templateimg = imresize(templateimg, [objsize(1)-1 tempsize(2)]);            
        elseif tempsize(1) < objsize(1) & tempsize(2) >= objsize(2)
            templateimg = imresize(templateimg, [tempsize(1) objsize(2)-1]);
        elseif tempsize(1) >= objsize(1) & tempsize(2) >= objsize(2)
            bCompare = 0;
        end
 
        if bCompare == 1
        
            dtemplateimg = double(templateimg);

            [x1, x2, y1, y2, score, nremoved] = TempMatchC(dtemplateimg, dobjectimg);
            
            %template size may have changed so get again
            tempsize = size(templateimg);
            templatearea = tempsize(1)*tempsize(2);
        
            pscore = score/templatearea;            

            if pscore < bestscore & x1 > 0 
                bestscore = pscore;
                bestindex = i;
            
                bestx1 = x1;
                bestx2 = x2;
                besty1 = y1;
                besty2 = y2;          
            end   
        end
end








                                    
function ncount = CreateTrainingSet(v, fin, fbackground, bsubtract, outfile, ntrainimgs, nfindobjs, maxbestscore)

load(fin)

objsf = find(objs == nfindobjs);

len = length(objsf);

fullimgs = {};
minstepsize = 10;

nstepsize = round(len/(ntrainimgs*10));
if nstepsize < minstepsize
    nstepsize = minstepsize;
end

ncount = 0;
i = 1;

nrows = 80;
ncols = 80;

while i < len & ncount < ntrainimgs
    index = objsf(i);
    vindex = index + nprocessstart;    
    
    if index < len
        img = read(v,vindex);        
        
        j = 1;
        
        while j <= nfindobjs & ncount < ntrainimgs
                            
            x1 = vx1s(index, j); 
            x2 = vx2s(index, j); 
            y1 = vy1s(index, j); 
            y2 = vy2s(index, j); 
            
            oimg = img(x1:x2, y1:y2, :);
                     
            badd = 1;
                
            if ncount > 0
            
                lenimgs = length(fullimgs);
                        
                im1 = imresize(oimg, [nrows ncols]);
                    
                k = 1;            
            
                while k <= lenimgs & badd == 1
                
                    im2 = cell2mat(fullimgs(k));                                                                               
                    im2 = imresize(im2, [nrows ncols]);                                          
                
                    [bestx1, bestx2, besty1, besty2, bestscore] = TempMatchSameSizeImages(im1, im2);
                           
                    if bestscore < maxbestscore
                        badd = 0;
                    end
                                            
                    k = k + 1;
                end            
            end 
                                 
            if badd == 1                                
                ncount = ncount + 1;
        
                rimgs(ncount) = {oimg(:,:,1)};
                gimgs(ncount) = {oimg(:,:,2)};
                bimgs(ncount) = {oimg(:,:,3)};            
        
                fullimgs(ncount) = {oimg};
            
                sz = size(oimg);
            
                xlengths(ncount) = sz(1);
                ylengths(ncount) = sz(2);
            end
           
            j = j + 1;
        end
                  
    end
    
    i = i + nstepsize;
end



save(outfile, 'rimgs', 'gimgs', 'bimgs', 'fullimgs', 'xlengths', 'ylengths');



function [bint, l, r, t, b, poverlap] = GetOverlap(l1, r1, t1, b1, l2, r2, t2, b2)


bint = 0;

l = 0;
r = 0;
t = 0;
b = 0;

poverlap = 0;

if (l1 <= r2 & l2<=r1 & t1<=b2 & t2<=t1) | (l2 <= r1 & l1<=r2 & t2<=b1 & t1<=t2)...
   | (l1 >= l2 & r1 <= r2 & t1 >= t2 & b1 <= b2) | (l2 >= l1 & r2 <= r1 & t2 >= t1 & b2 <= b1) 
    bint = 1;
end

if bint == 1
    l = max(l1, l2); 
    r = min(r1, r2);
    t = max(t1, t2);  
    b = min(b1, b2);
    
    area1 = (r1-l1)*(b1-t1);
    area2 = (r2-l2)*(b2-t2);
    
    area = (r-l)*(b-t);
    
    
    poverlap = min(area/area1, area/area2);
end
        
